#include "math/math.h"
#include "geometry/2d/vector.h"
#include "geometry/2d/geometry.h"
#include "delaunaytriangulation.h"

namespace fem
{
    //============================================================================================================================
    //============================================================================================================================
    DelaunayTriangulation::DelaunayTriangulation(const QStringList &f_ids) : m_CacheStep(0.02), m_BorderSegmentsBCInfo(f_ids.size())
    {
        for(int i = 0; i < f_ids.size(); ++i)
        {
            Q_ASSERT_X(!m_FunctionsIDTable.contains(f_ids[i]), "DelaunayTriangulation::DelaunayTriangulation", "function id isn't unique");
            m_FunctionsIDTable[f_ids[i]] = i;
        }

        m_Cache = new QList<Triangle*>*[60];
        for(int i = 0; i < 60; ++i)
            m_Cache[i] = new QList<Triangle*>[60];
    }

    //============================================================================================================================
    //============================================================================================================================
    void DelaunayTriangulation::addBorderSegments(int superstructure, const QList<int> &segments)
    {
        for(int i = 0; i < segments.size(); ++i)
        {
            if(!m_Superstructures[superstructure].Points.contains(m_BorderSegments[segments[i]].firstPoint()))
                m_Superstructures[superstructure].Points.append(m_BorderSegments[segments[i]].firstPoint());

            for(int j = 1; j < m_BorderSegments[segments[i]].numOfPoints() - 1; ++j)
                m_Superstructures[superstructure].Points.append(m_BorderSegments[segments[i]].point(j));

            if(!m_Superstructures[superstructure].Points.contains(m_BorderSegments[segments[i]].lastPoint()))
                m_Superstructures[superstructure].Points.append(m_BorderSegments[segments[i]].lastPoint());
        }
    }

    //================================================================================================================================
    //================================================================================================================================
    int DelaunayTriangulation::addBorderSegment(int begin, int end, int num_of_knots)
    {
        QList<int> point_indices;
        point_indices.append(begin);
        for(int i = 0; i < num_of_knots; i++)
        {
            double x = m_Points[begin].x() + (m_Points[end].x() - m_Points[begin].x())*(i + 1)/(num_of_knots + 1);
            double y = m_Points[begin].y() + (m_Points[end].y() - m_Points[begin].y())*(i + 1)/(num_of_knots + 1);
            m_Points.append(geometry::_2d::Point(x,y));
            point_indices.append(m_Points.size() - 1);
        }
        point_indices.append(end);

        m_BorderSegments.append(BorderSegment(point_indices));
        return m_BorderSegments.size() - 1;
    }

    //============================================================================================================================
    //============================================================================================================================
    void DelaunayTriangulation::addBorderSegmentBCInfo(const QString &f_id, int segment, int bc_type, int bc_formula)
    {
        Q_ASSERT_X(m_FunctionsIDTable.contains(f_id), "DelaunayTriangulation::addBorderSegmentExtraInfo", "function id is unknown");
        Q_ASSERT_X(bc_type >= 1 && bc_type <= 4, "DelaunayTriangulation::addBorderSegmentExtraInfo", "unknown bc type");
        Q_ASSERT_X(bc_type != 0 || bc_formula != -1, "DelaunayTriangulation::addBorderSegmentExtraInfo", "default bc_formula value for bc_type 1-3");

        BorderSegmentBCInfo bc_info = {segment, bc_type, bc_formula};
        m_BorderSegmentsBCInfo[m_FunctionsIDTable[f_id]].append(bc_info);
    }

    //============================================================================================================================
    //============================================================================================================================
    void DelaunayTriangulation::prepareTriangleSuperstructure(int superstructure)
    {
        int points[3] = {
            m_Superstructures[superstructure].Points[0],
            m_Superstructures[superstructure].Points[1],
            m_Superstructures[superstructure].Points[2]
        };
        geometry::_2d::Vector v1(m_Points[points[0]], m_Points[points[1]]);
        geometry::_2d::Vector v2(m_Points[points[0]], m_Points[points[2]]);

        // Если z-компонента векторного произведения больше нуля, то точки расположены против часовой стрелки. Нужно поменять их порядок.
        if(v1.x()*v2.y() - v1.y()*v2.x())
            math::permute(m_Superstructures[superstructure].Points[1], m_Superstructures[superstructure].Points[2]);

        Rib *ribs = m_Superstructures[superstructure].addRibs(6);
        addRib(points[0], points[1], ribs + 0, ribs + 1);
        addRib(points[1], points[2], ribs + 2, ribs + 3);
        addRib(points[2], points[0], ribs + 4, ribs + 5);
        Triangle *triangle = m_Superstructures[superstructure].addTriangle(ribs + 0, ribs + 2, ribs + 4);

        ribs[0].Neighbour = NULL;
        ribs[1].Neighbour = triangle;
        ribs[2].Neighbour = NULL;
        ribs[3].Neighbour = triangle;
        ribs[4].Neighbour = NULL;
        ribs[5].Neighbour = triangle;

        addToCache(triangle);
    }

    //============================================================================================================================
    //============================================================================================================================
    void DelaunayTriangulation::prepareQuadrangleSuperstructure(int superstructure)
    {
        int points[4] = {
            m_Superstructures[superstructure].Points[0],
            m_Superstructures[superstructure].Points[1],
            m_Superstructures[superstructure].Points[2],
            m_Superstructures[superstructure].Points[3]
        };
        geometry::_2d::Vector v1(m_Points[points[0]], m_Points[points[1]]);
        geometry::_2d::Vector v2(m_Points[points[0]], m_Points[points[2]]);

        // Если z-компонента векторного произведения больше нуля, то точки расположены против часовой стрелки. Нужно поменять их порядок.
        if(v1.x()*v2.y() - v1.y()*v2.x())
            math::permute(m_Superstructures[superstructure].Points[1], m_Superstructures[superstructure].Points[3]);

        Rib *ribs = m_Superstructures[superstructure].addRibs(10);
        addRib(points[0], points[1], ribs + 0, ribs + 1);
        addRib(points[0], points[2], ribs + 2, ribs + 3);
        addRib(points[1], points[3], ribs + 4, ribs + 5);
        addRib(points[2], points[3], ribs + 6, ribs + 7);
        addRib(points[0], points[3], ribs + 8, ribs + 9);
        Triangle *triangle_1 = m_Superstructures[superstructure].addTriangle(ribs + 2, ribs + 6, ribs + 9);
        Triangle *triangle_2 = m_Superstructures[superstructure].addTriangle(ribs + 1, ribs + 8, ribs + 5);

        ribs[0].Neighbour = triangle_2;
        ribs[1].Neighbour = NULL;
        ribs[2].Neighbour = NULL;
        ribs[3].Neighbour = triangle_1;
        ribs[4].Neighbour = triangle_2;
        ribs[5].Neighbour = NULL;
        ribs[6].Neighbour = NULL;
        ribs[7].Neighbour = triangle_1;
        ribs[8].Neighbour = triangle_1;
        ribs[9].Neighbour = triangle_2;

        addToCache(triangle_1);
        addToCache(triangle_2);
    }

    //============================================================================================================================
    //============================================================================================================================
    TriangleNet DelaunayTriangulation::triangulate()
    {
        for(m_CurStruct = 0; m_CurStruct < m_Superstructures.size(); ++m_CurStruct)
        {
            if(m_Superstructures[m_CurStruct].NumOfVertices == 3)
                prepareTriangleSuperstructure(m_CurStruct);
            else
                prepareQuadrangleSuperstructure(m_CurStruct);

            for(int j = m_Superstructures[m_CurStruct].NumOfVertices; j < m_Superstructures[m_CurStruct].Points.size(); ++j)
            {
                Rib *divided_rib;
                Triangle *triangle;
                triangle = findTriangle(m_Superstructures[m_CurStruct].Points[j], &divided_rib);

                if(divided_rib == NULL)
                {
                    Rib *ribs = m_Superstructures[m_CurStruct].addRibs(6);
                    addRib(triangle->Ribs[0]->Begin, m_Superstructures[m_CurStruct].Points[j], ribs + 0, ribs + 1);
                    addRib(triangle->Ribs[1]->Begin, m_Superstructures[m_CurStruct].Points[j], ribs + 2, ribs + 3);
                    addRib(triangle->Ribs[2]->Begin, m_Superstructures[m_CurStruct].Points[j], ribs + 4, ribs + 5);

                    Triangle *triangle_1 = m_Superstructures[m_CurStruct].addTriangle(triangle->Ribs[1], ribs + 4, ribs + 3);
                    Triangle *triangle_2 = m_Superstructures[m_CurStruct].addTriangle(triangle->Ribs[2], ribs + 0, ribs + 5);

                    ribs[0].Neighbour = triangle;
                    ribs[1].Neighbour = triangle_2;
                    ribs[2].Neighbour = triangle_1;
                    ribs[3].Neighbour = triangle;
                    ribs[4].Neighbour = triangle_2;
                    ribs[5].Neighbour = triangle_1;
                    triangle->Ribs[1]->Twin->Neighbour = triangle_1;
                    triangle->Ribs[2]->Twin->Neighbour = triangle_2;

                    m_Superstructures[m_CurStruct].setTriangle(triangle, triangle->Ribs[0], ribs + 2, ribs + 1);

                    testTriangle(triangle);
                    testTriangle(triangle_1);
                    testTriangle(triangle_2);
                }
                else
                {
                    if(divided_rib->Neighbour != NULL)
                    {
                        Rib *ribs = m_Superstructures[m_CurStruct].addRibs(8);
                        addRib(divided_rib->Begin,           m_Superstructures[m_CurStruct].Points[j], ribs + 0, ribs + 1);
                        addRib(divided_rib->End,             m_Superstructures[m_CurStruct].Points[j], ribs + 2, ribs + 3);
                        addRib(divided_rib->Next->End,       m_Superstructures[m_CurStruct].Points[j], ribs + 4, ribs + 5);
                        addRib(divided_rib->Twin->Next->End, m_Superstructures[m_CurStruct].Points[j], ribs + 6, ribs + 7);

                        Rib *temp_ribs[4];
                        temp_ribs[0] = divided_rib->Next;
                        temp_ribs[1] = divided_rib->Next->Next;
                        temp_ribs[2] = divided_rib->Twin->Next;
                        temp_ribs[3] = divided_rib->Twin->Next->Next;

                        Triangle *triangle_1 = m_Superstructures[m_CurStruct].addTriangle(temp_ribs[0], ribs + 4, ribs + 3);
                        Triangle *triangle_2 = m_Superstructures[m_CurStruct].addTriangle(temp_ribs[3], ribs + 2, ribs + 7);

                        ribs[0].Neighbour = divided_rib->Neighbour;
                        ribs[1].Neighbour = triangle;
                        ribs[2].Neighbour = triangle_1;
                        ribs[3].Neighbour = triangle_2;
                        ribs[4].Neighbour = triangle;
                        ribs[5].Neighbour = triangle_1;
                        ribs[6].Neighbour = triangle_2;
                        ribs[7].Neighbour = divided_rib->Neighbour;
                        temp_ribs[0]->Twin->Neighbour = triangle_1;
                        temp_ribs[3]->Twin->Neighbour = triangle_2;

                        m_Superstructures[m_CurStruct].setTriangle(triangle, temp_ribs[1], ribs + 0, ribs + 5);
                        m_Superstructures[m_CurStruct].setTriangle(divided_rib->Neighbour, temp_ribs[2], ribs + 6, ribs + 1);

                        testTriangle(triangle);
                        testTriangle(divided_rib->Neighbour);
                        testTriangle(triangle_1);
                        testTriangle(triangle_2);
                    }
                    else
                    {
                        Rib *ribs = m_Superstructures[m_CurStruct].addRibs(6);
                        addRib(divided_rib->Begin,     m_Superstructures[m_CurStruct].Points[j], ribs + 0, ribs + 1);
                        addRib(divided_rib->End,       m_Superstructures[m_CurStruct].Points[j], ribs + 2, ribs + 3);
                        addRib(divided_rib->Next->End, m_Superstructures[m_CurStruct].Points[j], ribs + 4, ribs + 5);

                        Rib *temp_ribs[2];
                        temp_ribs[0] = divided_rib->Next;
                        temp_ribs[1] = divided_rib->Next->Next;

                        Triangle *triangle_1 = m_Superstructures[m_CurStruct].addTriangle(temp_ribs[0], ribs + 4, ribs + 3);

                        ribs[0].Neighbour = NULL;
                        ribs[1].Neighbour = triangle;
                        ribs[2].Neighbour = triangle;
                        ribs[3].Neighbour = NULL;
                        ribs[4].Neighbour = triangle;
                        ribs[5].Neighbour = triangle_1;
                        temp_ribs[0]->Twin->Neighbour = triangle_1;

                        m_Superstructures[m_CurStruct].setTriangle(triangle, temp_ribs[1], ribs + 0, ribs + 5);

                        testTriangle(triangle);
                        testTriangle(triangle_1);
                    }
                }
            }
        }

        for(int i = 0; i < m_Superstructures.size(); ++i)
            for(int j = 0; j < m_Superstructures[i].Triangles.size(); ++j)
                m_TriangleElements.append(TriangleElement(m_Superstructures[i].Triangles[j]->Ribs[0]->Begin,
                                                          m_Superstructures[i].Triangles[j]->Ribs[1]->Begin,
                                                          m_Superstructures[i].Triangles[j]->Ribs[2]->Begin,
                                                          m_Superstructures[i].Material));
        return TriangleNet(m_TriangleElements, m_Points);
    }

    //============================================================================================================================
    //============================================================================================================================
    DelaunayTriangulation::Triangle* DelaunayTriangulation::findTriangle(int point, Rib **divided_rib)
    {
        QList<Triangle*> rect_content = m_Cache[(int)(m_Points[point].y()/m_CacheStep)][(int)(m_Points[point].x()/m_CacheStep)];

        for(int i = 0; i < rect_content.size(); i++)
        {
            geometry::_2d::Point A = m_Points[rect_content[i]->Ribs[0]->Begin];
            geometry::_2d::Point B = m_Points[rect_content[i]->Ribs[1]->Begin];
            geometry::_2d::Point C = m_Points[rect_content[i]->Ribs[2]->Begin];
            geometry::_2d::Point Z = m_Points[point];

            double s = geometry::_2d::squareTriangle(A, B, C);
            double s1 = geometry::_2d::squareTriangle(Z, B, C);
            double s2 = geometry::_2d::squareTriangle(A, Z, C);
            double s3 = geometry::_2d::squareTriangle(A, B, Z);

            if(qAbs(s - s1 - s2 - s3) > DBL_EPSILON)
                continue;
            else
            {
                if(qAbs(s1) < DBL_EPSILON) *divided_rib = rect_content[i]->Ribs[1];
                else if(qAbs(s2) < DBL_EPSILON) *divided_rib = rect_content[i]->Ribs[2];
                else if(qAbs(s3) < DBL_EPSILON) *divided_rib = rect_content[i]->Ribs[0];
                else *divided_rib = NULL;
                return rect_content[i];
            }
        }
        return NULL;
    }

    //============================================================================================================================
    //============================================================================================================================
    bool DelaunayTriangulation::testDelaunay(Rib* conj_rib)
    {
        if(conj_rib->Neighbour == NULL)
            return true;
        else
        {
            geometry::_2d::Point A = m_Points[conj_rib->Next->End];
            geometry::_2d::Point B = m_Points[conj_rib->Begin];
            geometry::_2d::Point C = m_Points[conj_rib->Twin->Next->End];
            geometry::_2d::Point D = m_Points[conj_rib->End];

            double a,b,c,d;
            a = qAbs((A.x() - B.x())*(A.y() - D.y()) - (A.x() - D.x())*(A.y() - B.y()));
            b = (C.x() - B.x())*(C.x() - D.x()) + (C.y() - B.y())*(C.y() - D.y());
            c = (A.x() - B.x())*(A.x() - D.x()) + (A.y() - B.y())*(A.y() - D.y());
            d = qAbs((C.x() - B.x())*(C.y() - D.y()) - (C.x() - D.x())*(C.y() - B.y()));

            if(a*b + c*d >= -DBL_EPSILON)
                return true;
            else
                return false;
        }
    }

    //============================================================================================================================
    //============================================================================================================================
    void DelaunayTriangulation::testTriangle(Triangle *triangle)
    {
        for(int i = 0; i < 3; ++i)
            if(!testDelaunay(triangle->Ribs[i]))
            {
                Rib *ribs = m_Superstructures[m_CurStruct].addRibs(2);
                addRib(triangle->Ribs[i]->Twin->Next->End, triangle->Ribs[i]->Next->End, ribs + 0, ribs + 1);

                Rib *r_ribs[4];
                r_ribs[0] = triangle->Ribs[i]->Next;
                r_ribs[1] = triangle->Ribs[i]->Next->Next;
                r_ribs[2] = triangle->Ribs[i]->Twin->Next;
                r_ribs[3] = triangle->Ribs[i]->Twin->Next->Next;

//                triangle->Ribs[i]->Next->Twin->Twin = ribs + 2;
//                triangle->Ribs[i]->Next->Next->Twin->Twin = ribs + 3;
//                triangle->Ribs[i]->Twin->Next->Twin->Twin = ribs + 4;
//                triangle->Ribs[i]->Twin->Next->Next->Twin->Twin = ribs + 5;

                Triangle *neighbour = triangle->Ribs[i]->Neighbour;

                ribs[0].Neighbour = neighbour;
                ribs[1].Neighbour = triangle;
                r_ribs[0]->Twin->Neighbour = neighbour;
                r_ribs[1]->Twin->Neighbour = triangle;
                r_ribs[2]->Twin->Neighbour = triangle;
                r_ribs[3]->Twin->Neighbour = neighbour;

                m_Superstructures[m_CurStruct].setTriangle(triangle, r_ribs[1], r_ribs[2], ribs + 0);
                m_Superstructures[m_CurStruct].setTriangle(neighbour, r_ribs[0], ribs + 1, r_ribs[3]);

                testTriangle(triangle);
                testTriangle(neighbour);
                return;
            }

        addToCache(triangle);
    }

    //============================================================================================================================
    //============================================================================================================================
    void DelaunayTriangulation::addToCache(Triangle *triangle)
    {
        int min_rect_index_v;
        int max_rect_index_v;
        int min_rect_index_h;
        int max_rect_index_h;

        int points[3] = {triangle->Ribs[0]->Begin, triangle->Ribs[1]->Begin, triangle->Ribs[2]->Begin};

        min_rect_index_v = qMin(m_Points[points[0]].x()/m_CacheStep, qMin(m_Points[points[1]].x()/m_CacheStep, m_Points[points[2]].x()/m_CacheStep));
        max_rect_index_v = qMax(m_Points[points[0]].x()/m_CacheStep, qMax(m_Points[points[1]].x()/m_CacheStep, m_Points[points[2]].x()/m_CacheStep));
        min_rect_index_h = qMin(m_Points[points[0]].y()/m_CacheStep, qMin(m_Points[points[1]].y()/m_CacheStep, m_Points[points[2]].y()/m_CacheStep));
        max_rect_index_h = qMax(m_Points[points[0]].y()/m_CacheStep, qMax(m_Points[points[1]].y()/m_CacheStep, m_Points[points[2]].y()/m_CacheStep));
        for(int j = min_rect_index_v; j <= max_rect_index_v; j++)
            for(int k = min_rect_index_h; k <= max_rect_index_h; k++)
                // Это очень медленная операция. Оказалось, что намного быстрее иногда проверять принадлежность точки к одному и тому же
                // треугольнику по нескольку раз, чем обеспечивать не повторяемость треугольников в кэше.
                //if(!m_Cache[k][j].contains(triangle))
                    m_Cache[k][j].prepend(triangle);
    }

    //============================================================================================================================
    //============================================================================================================================
    void DelaunayTriangulation::getDirichletBoundaryInfo(const QString &f_id, QList<int> &points, QList<int> &conditions) const
    {
        Q_ASSERT_X(m_FunctionsIDTable.contains(f_id), "DelaunayTriangulation::getDirichletBoundaryInfo", "function id is unknown");

        points.clear();
        conditions.clear();

        for(int i = 0, index = m_FunctionsIDTable[f_id]; i < m_BorderSegmentsBCInfo[index].size(); i++)
        {
            if(m_BorderSegmentsBCInfo[index][i].BCType == 1)
            {
                int segment_id = m_BorderSegmentsBCInfo[index][i].Segment;

                if(!points.contains(m_BorderSegments[segment_id].firstPoint()))
                {
                    points.append(m_BorderSegments[segment_id].firstPoint());
                    conditions.append(m_BorderSegmentsBCInfo[index][i].BCFormula);
                }

                for(int j = 1; j < m_BorderSegments[segment_id].numOfPoints() - 1; j++)
                {
                    points.append(m_BorderSegments[segment_id].point(j));
                    conditions.append(m_BorderSegmentsBCInfo[index][i].BCFormula);
                }

                if(!points.contains(m_BorderSegments[segment_id].lastPoint()))
                {
                    points.append(m_BorderSegments[segment_id].lastPoint());
                    conditions.append(m_BorderSegmentsBCInfo[index][i].BCFormula);
                }
            }
        }
    }
}
